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Abstract. We perform a density-matrix renormalization-group study of strongly interacting 
bosons on a three-leg ladder in the presence of a homogeneous flux. Focusing on one-third 
filling, we explore the phase diagram in dependence of the magnetic flux and the inter-leg 
tunneling strength. We find several phases including a Meissner phase, vortex liquids, a 
vortex lattice, as well as a staggered-current phase. Moreover, there are regions where the 
chiral current reverses its direction, both in the Meissner and in the staggered-current phase. 
While the reversal in the latter case can be ascribed to spontaneous breaking of translational 
invariance, in the first it stems from an effective flux increase in the rung direction. Interactions 
are a necessary ingredient to realize either type of chiral-current reversal. 
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1. Introduction 

Experimental progress with ultracold quantum gases has made feasible engineering the 
coupling between the different states of the atoms, in order to realize synthetic gauge 
fields [1,2]. The effective magnetic fields acting on the neutral atoms can be much larger than 
what is possible in solid-state systems. These advances bring the simulation of a wide range 
of Hamiltonians into reach that are important in condensed matter physics [3, 4, 5, 6, 7, 8]. 
Indeed, some of the most intriguing phenomena in condensed matter physics involve the 
presence of strong magnetic fields. For instance, topological states of matter are realized 
in quantum Hall systems [9, 10], which are insulating in the bulk, but bear conducting edge 
states. Remarkably, topological phase transitions were observed in experiments with cold 
atoms [11, 12, 13]. 

Recently, there has been a growing theoretical (see, e.g., work on fermions [14, 15, 16] 
and bosons [17, 18, 19, 20, 21, 22, 23]) and experimental [24, 25, 26] interest in quasi-one 
dimensional relatives of the square lattice [27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37], the iV-leg 
ladder systems with synthetic gauge fields. Besides using superlattices [24] to realize these 
geometries, a synthetic lattice dimension can be exploited where the sites on the rungs of the 
ladder correspond to different hyperfine states [38]. The latter has motivated several recent 
theoretical studies [39, 40, 41, 42, 43, 44]. Chiral edge currents on two- and three-leg ladders 
subjected to a homogeneous flux have been observed using both experimental approaches, 
for bosons [24, 26] and fermions [25]. The physics of interacting bosons on ladders has also 
been studied in the context of Josephson junction arrays [45, 46, 47, 48], albeit in the weakly- 
interacting regime. 

Such iV-leg ladder systems, similar to their counterparts in quantum magnetism [49], 
provide an elegant bridge between the physics in one and two dimensions [50]. Interacting 
bosons on two-leg ladders in homogeneous magnetic fields harbor very rich physics, featuring 
Meissner and vortex phases [51, 17, 18, 19, 52, 53, 21, 23], a biased-ladder phase that breaks 
the symmetry between the two legs [54, 55, 22], analogues of Laughlin states [56] as well as 
an interaction-driven reversal of the direction of the current due to spontaneous breaking of 
translational invariance in vortex lattices [22] . 

Here, we use density-matrix renormalization-group (DMRG) [57, 58, 59] simulations 
to explore the phase diagram of strongly interacting bosons on a three-leg ladder subjected 
to a homogeneous flux. We focus on a filling of one-third of a boson per site, which can 
easily be realized in experiments [26]. At this filling, the system is expected to be in a 
Mott-insulating state, based on the existence of a magnetization plateau at one-third of the 
saturation magnetization in the closely related three-leg spin-1/2 ladders [60, 61] and on work 
on the MI-SF transition in bosonic three-ladders at zero flux a zero flux [62]. We present the 
phase diagram for hard-core bosons (HCBs) in dependence of the inter-leg coupling strength 
and the magnitude of the magnetic flux. As a main result, we observe vortex phases that 
cannot be traced back to features in the single-particle dispersion. In the Meissner phase 
in this system, there is a fascinating reversal of the direction of the current, driven by the 
magnetic flux. We provide an explanation for this effect and we show that it persists at 
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Figure 1 . Schematic representation of the three-leg ladder model (1). (r,£) labels the site 
on the r-th rung and on the £-th leg. The tunneling strength along the legs(rungs) is ./(,/ ). 
We choose the gauge such that the phase is picked up on the rungs, resulting in a net flux of <b 
through each plaquette. The on-site interaction strength is given by U. 


intermediate interaction strengths. Finally, a staggered-current phase, generalizing the one 
found on the two-leg ladder [17, 18], is observed around 0 ~ ir, which also triggers a chiral- 
current reversal. 


2. Model and method 

We consider the Bose-Hubbard Hamiltonian on a three-leg ladder of length L (see figure 1): 

L 2 

K = -J ± J2Y1 ( e ~ ir * a t,l a r,l +1 + h.C.^j ( 1 ) 

r=l e= 1 

L -1 3 U L 3 

- J EE( a r,e a r+1,£ 3 " h.C^j + — EE", - 1 ), 

r=l 1=1 r=l 1=1 

where a\ t creates a boson on the r-th rung and the f-th leg of the system and n r y = a\ £ a r y 
is the on-site number operator. J and ./ are the tunneling matrix elements along the legs and 
rungs, respectively, U is the on-site interaction strength, and 0 is the magnetic flux through a 
plaquette of the ladder. 

We compute the ground state of (1) numerically with a single-site DMRG algorithm [63]. 
While exploiting the U( 1) symmetry of the Hamiltonian associated to particle number 
conservation, we keep up to 4000 DMRG states and the data in the figures of the main text are 
for L = 100, while we have considered L as large as L = 200 in the appendix. For U/ J < oo, 
we limit the number of bosons per site to six at maximum. For the determination of the Mott 
insulator (MI) to superfluid (SF) transitions, we use infinite-size DMRG [64] with up to 600 
DMRG states [see Appendix B]. 

The different phases are primarily characterized by their local current configurations. 
The associated operators are 

Jrj: tJ : 

Jr,l = iJ ± ( 




( 2 ) 

(3) 
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Figure 2. Phase diagram for HCBs with 1/3 particle per site on average, as a function 
of J L and (b. The dotted lines are upper bounds for the Mott-insulator (MI) to superfluid 
(SF) transition, as determined by an analysis of the Binder cumulant corresponding to an 
appropriate string-order parameter [see Appendix B]. The phase diagram features Meissner 
(M-MI/M-SF), vortex-liquid (V), vortex-lattice (VL) and staggered-current (SC-MI/SC-SF) 
phases. While the SC region is split into a Mott insulating and a superfluid phase, the V and 
VL phases are Mott insulating. In the Meissner phase, the color encodes the amplitude of the 
chiral current j c . The gray ’+’ symbols mark the parameter points analyzed with DMRG. 
Large black diamonds denote the parameter values for which the local-current patterns are 
presented in figure 3. The vortex-liquid region contains one- and two-component Luttinger 
liquids. The VL is expected to be surrounded by a (possibly very thin) vortex-liquid phase. 


We compute the chiral current from 

Jc = -Jr,3>- W 

r 

Throughout the paper, we fix the filling to 1/3 bosons per site. 

3. Phase diagram 

3.1. Summary of results for two-leg ladders 

Let us first briefly summarize some results for the corresponding two-leg ladder model in the 
low density regime (obtained by limiting the number of legs to two in figure 1), which will 
help to appreciate our results for the three-leg ladder to be presented in the following. The 
phase diagram for hard-core bosons at a filling of one boson per rung has been reported 
in [21, 23]: at fixed rung coupling J ± /J < 1.6, there is a transition from a phase with 
Meissner-like currents to a vortex phase when increasing the flux. Above a critical value 
J ± c ~ 1.6 J, the Meissner phase is stable at any flux. Furthermore, at certain commensurate 
vortex densities, vortex lattices are expected to form [51], which, however, do not seem to exist 
in the limit of hard-core bosons [21]. At smaller interaction strength, though, stable vortex 
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lattices at vortex densities of 1/2 and 1/3 have been identified in [22]. In those phases, the 
translation invariance associated with translations along the legs of the ladder is spontaneously 
broken and an enlarged unit cell forms. As a consequence, the effective flux sensed by the 
bosons is modified. This leads to a spontaneous reversal of the direction of the chiral current, 
under certain conditions on the flux and interaction strength [22] . 

3.2. Overview of the phase diagram for the three-leg ladder 

Figure 2 shows the ground-state phase diagram for HCBs [( a ) ri ) 2 = 0; U/J —> oo] of the 
three-leg ladder model equation (1) as a function of J ± e [J, 2J] and flux. Results are 
displayed for <f> e [0,7r] since all physical quantities are 27r-periodic, odd or even, functions 
of the flux [22]. We find the following phases: (i) a Meissner phase (M-MFM-SF), which 
shows a reversal of the current direction for large values of the flux, (ii) vortex-liquid phases 
(V), (iii) a vortex-lattice phase (VL), and (iv) staggered-current phases (SC-MI/SC-SF). The 
corresponding transitions are located by solid lines in Fig. 2, in which the gray ’+’ symbols 
indicate the parameter points that were analyzed numerically. For J ± /J > 1.6 or o > 0.37T, 
we observe a finite mass gap [see Appendix B] indicating that the system is in a Mott- 
insulating state. This is inherited from the J ± /J 1 limit, in which the ground state is the 
product of the local ground state on each rung, separated by a gap V2J, from all excited states 
(similar to the two-leg ladder with one boson per rung [65]). Upon lowering J ± at fixed o, the 
mass gap eventually closes. The computation of the Binder cumulant [66] of an appropriate 
string-order parameter for the Mott insulator yields an upper bound for the transition to a 
superfluid phase when lowering J x (see the dotted line in figure 2). This transition is, for 
o/x < 0.3, compatible with a Berezinskii-Kosterlitz-Thouless transition [see Appendix B] as 
expected from the theoretical work on magnetization plateaux in A'-leg spin ladders [60, 61] 
and numerical work at <f> = 0 [62]. The transition at 0.9 < f/n < 1, on the other hand, 
is compatible with a second-order phase transition [see Appendix B]. In the following, we 
refer to states with a Meissner-like current configurations using the label M, and to states with 
a staggered-current pattern using the label SC (irrespective of whether these are MI or SF 
phases). 

3.3. Meissner phases 

For J ± /J > 1.6, only a Meissner phase exists, characterized by currents occurring exclusively 
in the upper and the lower leg [see figures 3(a) and 3(e)], with a finite mass gap. The 
existence of such a Meissner-Mott insulator (M-MI) on the three-leg ladder has been predicted 
in reference [56]. For 4>/tt < 0.3, we find a transition from this Mott-insulator to a superfluid 
phase with Meissner currents (M-SF) by lowering the inter-chain coupling J ± [dashed line 
in figure 2], Intriguingly, the chiral current reverses its chirality from counterclockwise 
[figure 3(a)] to clockwise [figure 3(e)] for f/n > 0.75, meaning that the atoms flow in the 
direction opposite to the one favored by the effective magnetic field. As a consequence, at the 
point at which the reversal occurs, no current flows even though the bosons feel a very strong, 
non-staggered, magnetic flux. A typical example for the j c = j c {4>) curve in the Meissner 
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Figure 3. Current patterns for HCBs in (a) the Meissner phase (J±/J = 2, cI/tt = 0.2); 
(b-c) in the vortex-liquid phases ( J ± /J = 1.2, (f>/ir = 0.375 and 0.4); (d) in the VL phase 
( J ± /J = 1.2, <j)/n = 0.6); (e) in the Meissner phase, where the chiral current is reversed 
( J±/J = 1.2, <fi/ir = 0.9); and (f) in the staggered-current phase ( J ± /J = 1, <f> = 7r). The 
width of the arrows is proportional to the current strength and normalized to the strongest 
current in the region displayed for each parameter point. 


phase is shown in figure 4(a) for J ± /J = 1.8. The curve is smooth as there is no phase 
transition. 

3.4. Vortex phases 

Upon lowering J ± to J ± / J < 1.5, the Meissner phase is split by vortex phases at intermediate 
values of the flux. Typical current patterns in the vortex phases, which, for open boundary 
conditions, exhibit non-zero rung currents (jf e ), are presented in figures 3(b)-(d). We also 
plot the chiral current j c [see figure 4(a)] as a function of d for ./,/./ = 1.2, together with 
the average absolute value of the rung current j R = £=12 \ (jf t ) |/2 L [see figure 4(b)] . The 

transition from the Meissner into the vortex phases is characterized by j R becoming non-zero 
and by kinks in the j c = j c ((p) curves (see the example of J L /J = 1.2). The vortex phases can 
be further divided into vortex-liquid (V) phases, which are incommensurate, gapless phases 
[see, e.g., figure 3(b)-(c)] and a vortex lattice (VL), which is fully gapped and forms at a 
commensurate vortex density [51] [see, e.g., figure 3(d)]. 

3.4.1. Vortex lattice. The commensurability of the phases can be unveiled by studying the 
spatial patterns of the rung currents {jf l=l ). In general, their Fourier transforms bear two 
typical wavelengths qA and q B (see the inset of figure 5, we choose Iq^l < \qA\)- Motivated 
by this observation, we define two vortex densities p^ ,B = q,\ji(r ] /(2n) (with a the lattice 
spacing), which are shown in figure 5 for J ± / J = 1.2. In the Meissner and VL phases, both 
vortex densities are commensurate with p“ = 0 and = 1/2 (with a = A, B), respectively. 
These two phases are fully gapped. In the case of the VL, we ascribe this behavior to the 
incompressibility of the vortex pattern: it costs a finite energy to add a vortex to the system, 
even in the thermodynamic limit. The transitions to the vortex-liquid phases are continuous 
commensurate-incommensurate transitions [51]. We therefore expect the vortex-liquid region 
to surround the VL phase everywhere, even though the proximity of various phase transitions 
renders it very difficult to resolve numerically. 
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Figure 4. (a) Chiral current j c [cf. equation 4] as a function of the magnetic flux <j> for 
J±/J = 1.2 and 1.8 and HCBs. Inset: Comparison with perturbation theory (PT) [equation 5] 
for J ± /J = 5. (b) Average rung current j R = e=i 2 10^)1/27/ as a function of flux </> for 
J ± /J = 1.2 and 1.8. 


3.4.2. Vortex liquids. The incommensurate vortex-liquid region (V) encompasses phases in 
which both vortex densities are incommensurate [e.g., at (ft/it = 0.7 in figure 5] as well as 
phases where one mode is commensurate (at either p“ = 0 or p" = 1/2) and the other is 
not [e.g., at f/n = 0.375 and 0.5 in figure 5]. This is fully corroborated by the study of the 
von Neumann entropy, yielding either a central charge c = 2 or c = 1 [see Appendix A] 
in the vortex-liquid phases, corresponding to two- and one-component Luttinger liquids, 
respectively (we do not distinguish between the c = 1 and c = 2 vortex-liquid phase in the 
figures). In principle, the study of current-current correlations could also permit to distinguish 
the different vortex phases. However, the superposition of the behaviours of the different 
components renders the analysis of correlations less conclusive for finite-size systems. 

We stress that the emergence of vortex phases for bosons on the three-leg ladder is a 
many-body effect: the minimum of the single-particle dispersion is always at zero momentum, 
corresponding to a Meissner phase with counterclockwise chiral-current. The vortex phases 
are thus not inherited from finite-momentum global minima in the single-particle dispersion 
relation, unlike in two-leg ladders [1 1, 52, 21, 22]. 

3.4.3. Staggered-current phase. Around J ± ~ J and for </> < tt, a phase with staggered 
currents (SC) emerges [see figure 3(f)]. It is characterized by strong currents going around 
each plaquette. This phase breaks translation invariance along the legs of the ladder and has a 
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Figure 5. Vortex densities pfl B as a function of magnetic flux & for J ± /J = 1.2 and 
HCBs. They are extracted from the Fourier transform of (jfg =1 ), which show two dominant 
wavelengths in the vortex phase as illustrated in the inset for tp/ir = 0.45. 


commensurate vortex density p“ = 1/2, corresponding to a spatial periodicity of two lattice 
sites. It can be visualized as two VL states on the upper and lower two-leg ladder, shifted by 
one site with respect to each other. This phase shows a Mott-insulator to superfluid transition 
when lowering J ± / J [see Fig. 2 and Appendix B]. At (f> = it, the Hamiltonian (1) is real, and 
therefore, time-reversal invariant. The SC phases then spontaneously break this symmetry 
[see figure 3(f)]: reversing all currents does not leave the pattern invariant; they can thus 
be viewed as a generalization of the chiral Mott-insulator and superfluid phases realized on 
two-leg ladders [17, 18]. 

3.5. Chiral-current reversal in the Meissner phase 

Let us further investigate the region f/n > 0.75 where the chirality of the chiral currents is 
opposed to the one favored by the bare flux. In the strong-rung limit J ± J, the three-site 
problem is solved analytically and standard perturbation theory (PT) can be applied. To first 
order in J/J ± , the chiral current is then given by 



(5) 


which implies a reversal of the current at f = arccos(—2/3) ~ 0.737r, independently of 
J ± . Equation (5) agrees very well with the DMRG calculations for ./,/./ = 5 [see the inset of 
figure 4(a)], and the parameter point at which the reversal occurs barely moves for J ± / J > 1.6 
[see figure 2]. 

The flux dependence j c oc sin(0) is directly inherited from that of the two-leg ladder in 
the Meissner phase [21]. The two terms of (5) can therefore be interpreted as stemming from 
two contributions, which are sketched in figure 6. The first term (left panel) corresponds to 
chiral currents flowing on the two sub-ladders (formed by the middle leg and the upper or 
lower leg) resulting in a zero net-current on the middle leg, whereas the second term (right 
panel) corresponds to particles propagating only along the top and bottom leg of the three- 
leg ladder. This can be thought of as a Meissner phase with an effective flux f eS = 20. For 
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Figure 6. Schematic representation of the two contributions to equation (5) for n/2 < cf> < tt. 


7t/2 < 0 < 7T, the second contribution is negative, which can lead to a total chiral current with 
a reversed chirality. This reversal is thus associated with a doubling of the effective flux along 
the rung direction. It can already be captured in the minimal model of two plaquettes (L = 2), 
in which a qualitatively similar reversal occurs in the strongly interacting limit. An increase of 
the effective flux also underlies the chiral-current reversal in the case of two-leg ladders [22], 
yet there, the increase results from spontaneous breaking of translation symmetry. 

The J ± J limit also permits to understand another numerical observation: the density 
is imbalanced between the legs, in all phases. PT indeed predicts a density twice as large on 
the middle leg than on the outer ones. This difference decreases when lowering J ± /J. 

3.6. Finite interactions U/ J < oo 

In the single-particle case (U/ J = 0), the Meissner phase has counterclockwise chirality for 
any choice of parameters. Therefore, there must be a reversal of the chiral current as the 
interaction strength increases. This is indeed the case as shown in figure 7, displaying j c as 
a function of U, for f/n = 0.8 and 0.85 and J±/J = 1.6. The reversal from a negative to 
a positive current occurs at a finite interaction strength U whose value depends on 0. For 
0/0r = 0.85, the system enters into a staggered-current phase at intermediate values of U, 
as indicated in the figure. Translational invariance is spontaneously broken in the SC phase 
and the unit cell comprises four plaquettes [see figure 3(f)]. For 0 e [37r/4, tt], the effective 
flux is 0 e ff = 40 G [—7T, 0] modulo 27 t, and the current is reversed. This realizes another 
instance of the chiral-current reversal due to spontaneously enlarged unit cells first presented 
in reference [22] . The presence of the SC phases stabilizes the current reversal down to much 
smaller values of U/ J compared to parameters for which the SC phase is absent (compare 
the data for 0/V = 0.8 and 0.85 shown in figure 7). Note that the VL phase leads to the 
same unit-cell enlargement. However, as it occurs for 0 e [7 t/2, 37t/ 4] in the HCBs case (see 
figure 2), one has 0 e fr = 40 e [0,7r] modulo 27 t, and the chiral current is not reversed. 

4. Experimental realizations 

The most straightforward experimental approach would be to use a superlattice to split a two- 
dimensional lattice into three-leg ladders [67]. This would create an energy offset between 
the middle and the outer two legs, which could easily be compensated by another superlattice 
with a two-site periodicity. For systems with a synthetic lattice dimension [26, 25, 38], the 
interaction in the rung direction is not on-site as considered here, but depends on the total 
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Figure 7. Chiral current j c as a function of on-site interaction U for ,/, /,/ = 1.6 and 

0/tt = 0.8 and 0.85. 



Figure 8. Chiral current j c [cf. equation 4] as a function of the magnetic flux (t> for 
J ± /J = 0.2, 0.4 and 1, for long-range interactions along the rungs [cf. equation 6] and 
U/J = 10. 


density rfl = Y2t=\ n r j in all three rung sites. Thus, the appropriate form is 

U x— > 

Hmt = —^2 nKnl - 1). ( 6 ) 

r =1 

Our DMRG results for U/J = 10 show that this type of long-range interactions suppresses 
vortex phases. This is qualitatively consistent with the study of the two-leg ladder, in 
which this type of interaction also favors the Meissner phase over the vortex phase [19, 56]. 
The chiral-current reversal in the Meissner phase at a filling of one-third survives at small 
JjJ < 0.3 [see figure 8]. 
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5. Summary 


We presented the phase diagram of strongly-interacting bosons on a three-leg ladder subjected 
to a homogeneous magnetic flux. We identified several phases, including vortex-liquid phases, 
a vortex lattice and a Meissner phase. Moreover, there is a state with staggered currents, which 
leads to a reversal of the chiral current due to the spontaneous increase of the unit cell, similar 
to the situation discussed in [22]. Fascinatingly, the Meissner phase also shows a chiral- 
current reversal when increasing the strength of the magnetic flux per plaquette, yet in the 
Meissner phase, translational invariance is not broken. The analysis of the strong-rung and 
large-interaction limit indicates that a doubled flux is experienced along the rung direction. 
Therefore, the chiral-current reversal in the Meissner phase is qualitatively different. We argue 
that interactions are a necessary ingredient to obtain these behaviors. 
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Appendix A. Entanglement entropy 

The entanglement entropy S v n is defined as the von Neumann entropy corresponding to a 
bipartition of the wavefunction into two subsystems A and B: 


S vN = - Tr p A In p A , 


(A.l) 


where p,\ = Tr b U) (U is the reduced density matrix of subsystem A. In conformal field 
theory, for a one-dimensional system with open boundary conditions and total length L, the 
von Neumann entropy of a subsystem of length r is given by [ 68 , 69] 



Q 

S V N(r) = 77 log 
6 


(A.2) 


where c is the central charge of the conformal field theory and g is a non-universal constant. 
The central charge measures the number of gapless modes in the system. 

In figure A1 we show the entanglement entropy S v n(t) for straight cuts through the 
three-leg ladder between rungs r and r + 1. By fitting (A.2) to the data we obtain estimates 
for the central charge. Oscillations of the entropy in the vortex phases make the fit very 
sensitive to the fitting region used. This applies, in particular, to the vortex-liquid phases. 
Nevertheless, the behavior of the central charge provides further insights into the nature of the 
commensurate-incommensurate transition. 
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Figure Al. The entanglement entropy S v jv as a function of the relative length of the 
subsystem r/L for L = 100. Symbols show DMRG results for (a) J ± / J = 2.0, (j)/n = 0.2 in 
the Meissner phase; (b) J ± /J = 1.2, (j)/i r = 0.375 in the vortex-liquid phase; (c) J ± /J = 1.2, 
(f >/it — 0.6 in the vortex-lattice phase; (d) J ± / J = 1.2, (fi/n = 0.7 in the vortex-liquid phase; 
(e) Jl/J = 1.2, (b/r = 0.9 in the Meissner phase with currents with reversed chirality; (f) 
J ± /J = 1.0, (f>/i r = 1.0 in the staggered-current phase. The red line is a fit of (A.2) to the 
DMRG data. The insets show the mass gap A Em as a function of inverse system length which 
is extrapolated to the infinite length limit with a linear fit. 


In the Meissner phase [figures A 1(a) and (e)] the values of the central charge are very 
close to zero in agreement with an entanglement entropy independent of the subsystem size, 
which corresponds to the area law of a fully gapped phase [70]. 

The vortex phases can be divided into vortex-liquid and vortex-lattice phases with 
incommensurate and commensurate vortex densities, respectively. As illustrated in figure 5 
of the main text there are two vortex densities that undergo a transition from the Meissner to a 
vortex-liquid and vortex-lattice phase with increasing flux f/iT and at low enough J ± / J. This 
picture is consistent with the behavior of the central charge. More precisely, at J x /J = 1.2 
and f/n = 0.375, we find one incommensurate vortex density in the Fourier transform of the 
rung currents and a central charge of c ~ 1 [figure Al(b)]. At flux of" = 0.7, there are two 
incommensurate vortex densities. The central charge is c ~ 2 in agreement with the presence 
of two gapless modes [figure Al(d)]. At f/n = 0.6, both vortex densities are commensurate 
and we find an entanglement entropy independent of subsystem size (up to oscillations due to 
the increased unit cell), or c = 0 in terms of the central charge [figure Al(c)]. 

This suggests that the modes, which are measured by the central charge, can be 
related to the dominating vortex densities in the Fourier transform of the rung currents. 
An incommensurate vortex density would be associated with a gapless mode while a 
commensurate vortex density would correspond to a gapped mode. 

In the staggered-current phase the vortex densities are commensurate, while the fit gives 
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a central charge c ~ 1. But the scaling of the entanglement spectrum indicates that it is two 
independent c = 0.5 modes, and also the scaling of the correlation function gives the scaling 
dimension A ~ 1/8 (not shown), consistent with c = 0.5. This indicates that the mass gap 
is closed, which is compatible with the analysis of the mass gap and the Binder cumulant 
presented in the next section. 

Appendix B. Mott insulator to superfluid transition 

The mass gap A E M is defined as 


A E m = i [E(N - 1) + E(N + 1)] - E(N), 


(B.l) 


where E(N ) is the total ground-state energy of the system with N particles. 

We generally observe a finite mass gap that becomes smaller when decreasing the inter¬ 
leg coupling ,/j. In the J ± = 0 limit, the mass gap is known to be closed. As already 
mentioned in the main text, in the f/ir = 0 limit, the transition to the superfluid phase at 
finite J ± /J can be expected to be of the Berezinskii-Kosterlitz-Thouless type [60, 61], and, 
therefore, is very hard to pinpoint numerically. 

In the insets of figure A1 we show the mass gap as a function of the inverse system size 
1 /L extrapolated with a linear fit to the thermodynamic limit. We find a finite mass gap in 
the Meissner phase for J/J x = 2, f/ir = 0.2 and J/J ± = 1.2, f/ir = 0.9 [figures Al(a) 
and (e)]. The mass gap closes in the Meissner phase when lowering .//./, . The vortex-liquid 
phases [figures A 1(b) and (d)] and the vortex-lattice phase [figure A 1(c)] have finite mass gaps 
with weak system-size dependence. The finite mass gap in the vortex-liquid phases confirms 
that the charge mode does not contribute to the values of the central charges in figures Al(b) 
and (d) but they are fully accounted for by the two other modes on the three-leg ladder that 
undergo the commensurate-incommensurate transition. 

The analysis of the central charge is more difficult in the regions where the mass gap can 
not be resolved (or is actually zero) with the system sizes studied in this work. This is the case 
in the Meissner phase for J ± /J < 1.6 and 0/7r < 0.3 and in the staggered-current phase [see 
the inset of figure A 1(f)]. 

Calculating energy gaps is a quite difficult calculation when the gap is small, so where 
an order parameter exists it is generally better to use it instead to detect the phase boundary of 
the Mott-insulating region. A suitable order parameter for a Mott phase is the non-local string 
order parameter [71], given by 



(B.2) 


where n\ = ^ =1 n rt g is the total particle density on rung r. 
This can be expressed in terms of a ‘kink’ operator, 


P(i) = IW-1) 1 -* , 

from which O p takes the form of a correlation function, 
°l = .. iif 11 <P(*)P(J')> 7 


(B.3) 


(B.4) 
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Figure Bl. Binder cumulant B rn for various number of states m, for (a) <p/n = 0.2 and (b) 
<j>/i r = 0.95, to analyze the Mott insulator to superfluid transition. In (a) the crossing point 
moves towards lower J ± / J as m increases, compatible with a Berezinskii-Kosterlitz-Thouless 
transition at J x /J < 1.25. In (b) the crossing point does not move significatively when m 
increases, compatible with a second-order transition at J ± /J ~ 1.035. The average rung 
current j R , as defined in the main text, is also displayed. It shows that the Meissner to SC 
transition occurs at J ± /J ~ 1.135. 


or, we can construct an extensive order parameter P = ffpp ('/■)■ We cannot evaluate ( P) 
directly, since the sign is indeterminate, however, we can calculate ( P 2 ), which is simply 
related to O' 2 , as 


0 2 p = ( P 2 )/L 2 , (B.5) 

where L is the length of the lattice. The great advantage of this construction is that the P 
operator is easy to construct for matrix-product-state (MPS) algorithms using the triangular 
matrix-product-operator (MPO) formulation [72, 64], and the expectation value of such 
triangular MPO’s can be evaluated directly in the thermodynamic limit [73]. This allows 
us to calculate the Binder cumulant [66] of P, 


B 


(P*) 

3(P 2 ) 2 ' 


(B.6) 


For an infinite matrix-product-state (iMPS), the expectation value of the n-th power of an 
operator P n is obtained as a degree n polynomial of the lattice size L, which is exact in 
the asymptotic large-L limit. Hence the quantity O 2 = (P 2 ) / L 2 is evaluated directly as the 
coefficient of the degree 2 component of (P 2 ). 

If we evaluate the Binder cumulant directly for an iMPS in the large L limit, then it 
simply probes where O 2 is zero or non-zero, so does not give any additional information. 
Since the coefficient of the n-th degree term in ( P n ) is simply O p , it follows that if O 2 is 
non-zero then the Binder cumulant in the large L limit is identically equal to 2/3. Similarly, 
if O 2 = 0, then the cumulant expansion of (P 4 ) reduces to 3(P 2 ) 2 and the Binder cumulant 
is identically zero. Hence we have a step function located at the point where O 2 becomes 
non-zero, which is what one would expect if one takes the large-L limit of a finite-size Binder 
cumulant. 





















Strongly interacting bosons on a three-leg ladder in the presence of a homogeneous flux 15 

However, for an iMPS the correlation length is always finite, so one might expect that 
finite-entanglement-scaling with the basis size m can be used instead of finite-size scaling. 
This is indeed the case, if we evaluate the polynomial ( P n ) at L — bf where £ is the 
correlation length, and b is some scaling factor. This is equivalent to calculating the order 
parameter over a finite section of size L = b£ of the infinite lattice. 

Since £ depends on the number of states m, we can plot a family of curves of the Binder 
cumulant B m for different values of m. For a finite system, the value of the Binder cumulant 
at a second-order critical point is independent of L (up to higher order corrections). For an 
iMPS, the analogous result is that the critical value of the Binder cumulant is independent of 
m. Note that for finite size systems the actual value of the Binder cumulant at the critical 
point is not universal and depends on the boundary conditions [74]. For an iMPS, the value of 
the m-dependent Binder cumulant depends on the chosen scale factor b. For the calculations 
here, we use b = 1. 

As an example, we show in figure Bl(a) the Binder cumulant as a function of ./ for 
(f )/tt = 0.2. The data is compatible with a Berezinskii-Kosterlitz-Thouless transition from the 
Meissner Mott-insulating region into the Meissner superfluid region as J ± /J decreases and 
permits to give an upper bound for the critical point, which is reported in figure 2 in the main 
text. For f/n = 0.95 [see figure Bl(a)], a Mott-insulator to superfluid transition also occurs 
in the SC phase and it is compatible with a second-order critical point at J ± /J — 1.035. 
The Meissner to SC phase transition, which can be characterized by a rise of the average 
rung-current j R , occurs at a distinct point J x /J ~ 1.135. 
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